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We generalize the concept of quantum phase transitions, which is conventionally defined for a 
ground state and usually applied in the thermodynamic limit, to one for metastable states in fi- 
nite size systems. In particular, we treat the one-dimensional Bose gas on a ring in the presence 
of both interactions and rotation. To support our study, we bring to bear mean-field theory, i.e., 
the nonlinear Schrodinger equation, and linear perturbation or Bogoliubov-de Gennes theory. Both 
methods give a consistent result in the weakly interacting regime: there exist two topologically dis- 
tinct quantum phases. The first is the typical picture of superfluidity in a Bose-Einstein condensate 
on a ring: average angular momentum is quantized and the superflow is uniform. The second is 
new: one or more dark solitons appear as stationary states, breaking the symmetry, the average 
angular momentum becomes a continuous quantity. The phase of the condensate can therefore be 
continuously wound and unwound. 

PACS numbers: 03.75.Hh, 03.75.Lm 
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I. INTRODUCTION 

One of the main interests of a many-body quantum sys- 
tem is to identify its phase diagram and quantum phase 
transitions (QPTs). The concept of QPTs Q is usually 
formulated in the ground state in the thermodynamic 
limit. Some criteria for classification of a QPT are the 
order of transition and the universality class. Ideally, 
these can be determined from experimental observables 
such as energy, susceptibility, and other bulk measures, 
and their derivatives with respect to the system param- 
eter that triggers the transition; more recently, certain 
forms of entanglement have been proposed as alternate 
measures. A suitably defined order parameter charac- 
terizes the symmetry properties of a QPT. A classical 
phase transition is driven by thermal fluctuations while 
a QPT is driven by quantum fluctuations, and can be 
most unambiguously understood at zero temperature. 

In this article we show that the conventional ways of 
defining QPTs in a Bose-Einstein condensates (BECs) 
can also be applied to phase transitions in metastable 
states in finite systems. The issue of QPTs in ex- 
cited states has been addressed previously in nuclear 
physics [l, Q , where nuclei are generically finite and sig- 
natures of QPTs can be identified in certain excitation 
spectra. We choose a one-dimensional integrable sys- 
tem of bosons as a concrete example. Studies of one- 
dimensional systems relevant to our chosen model have 
a long history, including exactly solvable quantum sys- 
tems [4j46| , decay of persistent current |7Hll| , and classi- 
cal solitons [l^ . 

The most important fact concerning BECs in the ther- 
modynamic limit is Hohenberg's theorem |13j . which 
states that in an interacting and infinite homogeneous 
system condensation does not occur in less than three 



spatial dimensions at finite temperature and in one di- 
mension at absolute zero. This is specifically proven in 
the Lieb-Liniger model in the thermodynamic limit p^ . 
However, this is not true in finite-size systems or spa- 
tially confined systems motional ground-state con- 
densation in two and three dimensions is shown to be 
possible [iHl- Furthermore it is also proved via the 
Bogoliubov inequality that off-diagonal long-range order 
is allowed to exist not only in the ground state, but also 
in the general excited states. The exact theoretical treat- 
ment on low-dimensional Bose gases in finite-temperature 
is described for example in Ref. [l^ . 

In areas of study such as quantized vortices in super- 
fluid systems, metastability plays a key role in the study 
of quantum dynamics and transport properties. This is 
particularly the case for well-insulated systems, such as 
ultracold quantum gases, where there is only negligible 
exchange of energy and particles with the environment. 
In particular, in metastable states of matter waves, such 
as soliton trains [l^, [l^l , the effects of dissipation can be 
suppressed and a metastable condensate is observable. 
The main focus of past studies of QPTs has been on 
ground states and the thermodynamic limit, but ultra- 
cold quantum gases necessitate a reexamination of the 
role of excited states. Moreover, experiments in these 
systems have precise control over the effective dimension- 
ality, so that one and two dimensions can be studied for 
a wide range of interactions [2l| - [23| . both repulsive and 
attractive. 

In our previous analysis [l^l we showed that the aver- 
age angular momentum of weakly-repulsive bosons in a 
one-dimensional ring undergoes a continuous change in 
its metastable excited states as a function of interaction 
and rotation. This phenomenon is intuitively understood 
in terms of bifurcation of stationary excited-state energy 
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branches of the plane- wave state propagating on the ring, 
and of localized soliton trains within the mean-field the- 
ory. Each excited state has a denumerably infinite num- 
ber of bifurcations from the plane wave to a state con- 
taining 1,2,... dark solitons; each such bifurcation cor- 
responds to a QPT. Moreover, for attractive interactions 
there is a set of QPTs even for the ground state [25|. For 
completeness, we compare both repulsive and attractive 
interactions in our present mean-field study, although our 
main focus is on the repulsive case. In our previous anal- 
ysis [23 | we did not discuss in which parameter regime 
and in which energy regime the above phenomenon is 
observed. The aim of the present work is to reveals prop- 
erties of metastable QPTs within the mean-field theory 
and investigate the linear stability of the mean-field solu- 
tions in order to show their metastability. We specifically 
obtain the phase boundary for general soliton solutions, 
which are characterized by the number of density notches 
and phase winding number, in the experimentally vari- 
able parameters, i.e., the strength of repulsive interaction 
and the frequency of the rotating drive. The energy dia- 
gram between the different phases is also shown. In order 
to show that these excited-state solutions do not undergo 
any dynamical instability, and thereby that the solutions 
are metastable, we investigate the linear stability of all 
the solutions of the mean-field theory. 

Our presentation is organized as follows. In Sec. |TT] 
we introduce the Hamiltonian of our system, which cor- 
responds to the Lieb-Liniger model in a rotating frame, 
and describe its basic properties. In Sec. IIIII we derive 
all stationary states within the mean-field theory in the 
weakly interacting regime, and predict the transition be- 
tween the persistent-current state and localized soliton 
trains. In Sec. IIVI we discuss the linear stability of 
these solutions in order to determine their metastability 
against small perturbations. We summarize our results 
and discuss experimental possibilities to realize our ideas 
in Sec. El 



II. THE MODEL 

A. Lieb-Liniger Hamiltonian in a Rotating Frame 

Let us consider the Hamiltonian for one-dimensional 
bosons with a contact interaction, known as the Lieb- 
Liniger model Q , in position representation. 
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(1) 



where N is the number of bosonic atoms and giD is 
the effective strength of s-wave interatomic interaction in 
one dimension (ID) [2^. We impose a periodic bound- 
ary condition by assuming a ring-shaped waveguide or 
toroidal trap [27h33| of radius R. The coordinates (az- 
imuthal angles) of bosons are specified by a set of vari- 
ables {9} = {01, 92,... On} where 9j £ [0,2tt) for ah 



j € {!,..., iV}. The length and energy units are R 
and h^/{2mR^), respectively, and m is the mass of a 
boson. The coupling strength giD is measured in unit 
of h^/{2mR) and hence dimensionless. The Hamiltonian 
we address in this paper is the Lieb-Liniger Hamiltonian 
in the rotating frame of reference. 



where 
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(3) 



is the angular-momentum operator measured in units of 
h, the trap rotates at angular frequency f2, and the last 
term in is a constant energy associated with rigid- 
body rotation which is added to the Hamiltonian to make 
the system translationally invariant; however it does not 
change the results of this paper. 

As we will show later, states having angular momen- 
tum equals to an integral multiple of TV are always the 
stationary solutions (uniform supcrflow) for repulsive in- 
teractions of the Hamiltonian ([T} . In terms of the nonlin- 
ear Schrodinger equation approach, there exists another 
solution having a density modulation, so called "soli- 
ton" solutions. In the rest frame, the soliton solutions 
have node(s) where the density becomes zero, and such 
a density-modulated state with zero-density node(s) is 
called "black" or "dark" soliton. The uniform-density 
state and black-soliton state have an finite energy gap 
and they cannot cross over each other at finite strength 
of repulsive interaction. By the addition of the rotating- 
drive term {^2QL), in contrast, a density- modulated 
state is allowed to have nonzero density notch(es), which 
is called "gray" soliton. This fact makes continuous en- 
ergy change between distinct topological states, namely, 
uniform superflow and soliton states possible, because the 
phase modulation of the soliton state can range from in- 
finitesimal (almost constant phase) to maximum (corre- 
sponding to zero density notches). This is the main find- 
ing of our work, and we analytically derive this crossover 
in this paper. The terms gray and black are derived from 
density contrast ima gini ng in BECs and intensity imag- 
ing in optical fibers j34| . We use the term dark soliton 
train to refer to more than one equally spaced soliton, 
with gray or black indicating whether or not the solitons 
have nodes. 



B. Periodicity and Umklapp Processes 

All the physical quantities of Hamiltonian ([2]) are pe- 
riodic with respect to Q with period 1 in our units. In 
order to show this, we consider a Schrodinger equation 
as follows 133: 



Him{{9})^£in)^m). 



(4) 
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The many-body wave function ^{{9}) satisfies the single- 
valuedness boundary condition 



..,9n)^'^{9i,...,9j+2tt,...,9n) (5) 



for all atomic positions 9j. Substituting the transforma- 
tion 



*o({n) =cxp 



(6) 



into the Schrodinger equation ([4]) we can eliminate the 
fi-depcndent terms from the equation to yield 



(7) 



The boundary condition of the wave function is then 
modified to be 



,...,9j 
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Noting the fact that the new Schrodinger equation ([7]) 
does not depend on the angular frequency fl, and that 
the boundary condition ^ is periodic with respect to il 
with period of 1, we observe that the eigenvalue £o{^) 
must be periodic with the same period. Equation ^ 
states that once the eigensolutions of the Hamiltonian in 
the rest frame Hq are found, the wavefunction ^({0}) of 
H is obtained via the inverse transformation of Eq. ^ , 
and the energy is given by 



£ ^£o- 2n{L) +n^N , 



(9) 



where (L) = (4'o|I/|^'o) is the total angular momentum; 
this is a conserved quantity and thus a good quantum 
number. The periodicity with respect to fl enables us 
to understand our system in terms of the Bloch energy 
structure in a solid, although there is not any direct phys- 
ical relevance to band theory. In a solid, the energy dis- 
persion is periodic with respect to quasimomentum and 
one can thus define a reduced Brillouin zone in which all 
information of a system is included. In our periodic ring, 
the energy is periodic with respect to and one can re- 
strict the analysis within a primitive Brillouin zone [0, 1). 
Our restriction of € [0, 1) for the rest of this study is 
therefore completely general. 

The periodicity with respect to 51 naturally leads us to 
define the Umklapp process of the total momentum L' = 
L+JN with J G Z an arbitrary integer 0. The Umklapp 
excitations are translations of the center of mass on the 
ring. Namely, a state with angular momentum L has 
an infinite number of counterparts at angular momenta 
separated by integer multiples of the number of particles, 
each counterpart having identical properties. 



C. Conserved Quantities 

The reduced single-particle density matrix of an N- 
body wave function "^{{9}) is given by 



Pi I 



■2tt 



d92...d9N-^*{9',92,...,9N) 
x^f{9,92,...,9N)., 



The spectral decomposition of pi takes the form 



(10) 



(11) 



where ipj are the eigenvectors in the many-body Hilbert 
space and Xj are the associated eigenvalues. If there is 
one and only one dominant eigenvalue of pi , there exists 
off-diagonal long-range order (a Bose-Einstein conden- 
sate, or BEG), and the corresponding eigenfunction. 



(12) 



is regarded as an effective one-body quantity called the 
condensate wavefunction. Although Bose condensation 
does not occur in the thermodynamic limit in ID, here 
we consider finite systems for which TV does not tend to 
infinity, and the concept of BEG in ID is therefore valid. 

The superfluid velocity is defined from the phase of the 
condensate wavefunction in Eq. (E 



(13) 



v,{9) ^ -^^AiiO). 
m o9 



An integral over the superfluid velocity along a closed 
path is called the circulation, 



C= d9vs{9), 
Jo 



which is quantized as 



m m 



(14) 



(15) 



where m is the mass of the constituent particles in the 
condensate. 

Number, energy, and angular momentum are also con- 
served in the full quantum theory according to the usual 
relations. In addition, in the mean field theory there 
is a denumerably infinite set of conserved quantities as 
described by Zakharov [s^. The original many-body 
Hamiltonian ^ in one dimension is also integrable, e.g. 
via the Bethe ansatz. This will be addressed in a later 



work 37 1 . 



III. MEAN-FIELD SEMI-CLASSICAL THEORY 

A. Stationary Solutions 

When the contact atomic interaction is very weak, 
specifically, giriN < 0(1), the bosons form a conden- 
sate, whose static and dynamical properties are described 
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by the nonlinear Schrodinger, or Gross-Pitaevskii (GP) 
equation. It is very useful to study the stationary solu- 
tions of the GP equation to develop an intuitive knowl- 
edge of the condensate properties. 

For convenience we introduce a dimensionlcss param- 
eter 



7 



gipN 

2tt 



(16) 



that represents the ratio of the mean-field interaction en- 
ergy to the kinetic energy (note that our 7 is not the 
Tonks-Girardeau parameter). The GP equation under a 
rotating drive takes the form 



d 



where ip{9) is the condensate wavefunction from Eq. ([T2|) . 
Equation (fT7|) has two kinds of solutions (s^l: plane wave, 



iJ0 



and soliton train, 



(18) 



(19) 



where J € Z is the phase-winding number, and j € |Z| is 
the number of density notches in the soliton train. The 
derivation of the soliton solutions is summarized in Ap- 
pendix A. The amplitude of the soliton-train solution for 
repulsive interactions, giD > 0, is given by 



yp^ = AAJl + r/dn2 



jKie 



(20) 



where dn(u, k) is the Jacobi dn function with elliptic 
modulus k G [0,1], and K{k) and E{k) are elliptic in- 
tegrals of the first and second kind, respectively. The 
parameter is an arbitrary coordinate in the interval 
< ^0 < 27r, indicating that the soliton solution is a 
spontaneous-broken-symmetry state. 

The normalization constant Af is determined so that 
the order parameter is normalized to be unity on the ring: 



K 



2t:{K + f]E) 



(21) 



The depth 77 of density notches is obtained by substitu- 
tion of the amplitude into the GP equation as 



e[-l,0]. 



We define the functions 



/ = 7r27-2(jX)2 
g ^ f + 2{jKf, 
h = f + 2{jkK)\ 



2]'KE, 



(22) 



(23) 
(24) 
(25) 



for simplicity. Then the integration of the imaginary part 
of the GP equation gives an analytical expression for the 
phase prefactor: 



^g(0) = fir 




nU; 



jK{0 - 60) 



,fc , (26) 



where n(^, u, k) is an elliptic integral of the third kind 
with an amplitude parameter 



f 



and a sign function 



S ^ sign(C! - J) = ^ : J ^^j^ 



(27) 



(28) 



The soliton-train solution ([20)) has two limiting behav- 
iors. First, in the limit — > 0, both the amplitude and 
phase approach the plane-wave solution with the same 
phase winding J. Second, in the limit 77 — > — 1 (equiva- 
lent to / ^ 0), the wave function is found to approach 
the Jacobi sn function, which corresponds to a black soli- 
ton train with tt phase jumps and density notches which 
form nodes [s^- In this limit 77 — > — 1, the condensate 
wave function, chemical potential, and energy are given 
by 



(sn) 



^.(sn) 
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jXy (1 + k^)E-{l + 2P)K 
3iE-K) 



(30) 



(31) 
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FIG. 1: Amplitude |V'j'j'(^)l (upper panels) and phase pro- 
files ifijj (9) (lower panels) of the solitonic condensate wave- 
function Eq. (|19|) for various numbers of density notches 
j. Left: (i,7,n) = (1,0.7,0.45); middle: (ji',7,n) = 
(2,0.7,0.55); right: (j, 7, n) = (3, 0.7, 0.45). For the weakly- 
interacting regime the phase winding number J is determined 
by a fixed set of parameters {j, 7, Q,). 
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FIG. 2: (Color online) Parameter space where soliton solu- 
tions can exist. Filled regions present the regimes where the 
stationary (a) j = 1, (b) j = 2, and (c) J = 3 soliton solu- 
tions coexist with plane-wave solutions. Boundaries are given 
by the parabolas of Eq. (|56|) with a phase-winding number 
J, as indicated by an integer value next to the corresponding 
parabolic curve. The path C in the panel (c) denotes a typi- 
cal path to observe the continuous change in the topology of 
the order parameter (see, Sees. HITCl and UlI D|) . (d) Elliptic 
parameter for a single soliton on the ring {j = 1) as a func- 
tion of rotational drive Q £ [0, 0.5] and mean field strength 
7 G [—1.5, 1.5]. For Q = 0.5 the soliton solution is given by 
the Jacobi sn and cn functions for repulsive and attractive 
interactions, respectively. For SI = the soliton solution is 
given by the Jacobi dn function for attractive interactions. 



From analysis of the equation / = that determines the 
elliptic modulus, can be expanded near the critical 
point as 



4 10 2 
—7 - —7 
J J 



Oil' 



(32) 



In between these limits (— 1 < < 0) we say that Eq. (120)) 
describes gray solitons. In Fig. [T]we show typical ampli- 
tude and phase profiles of gray soliton trains for j ~ 1,2, 
and 3. 



B. Phase Diagram 

The value of the elliptic modulus k has so far been left 
undetermined. To practically obtain the physical quan- 
tities such as energy, the density profile and so on, one 
needs a concrete value of k. The elliptic modulus is deter- 
mined from the requirement that the phase of the order 
parameter satisfies the single- valuedness condition: 



(Pj^j{0 + 2tt) = ip.j,j{0) + 2TrJ, 



(33) 



boundary condition for the soliton. Equation (|33p can be 
read as 

27r\n~ J\ 




^+j^[l-Ao(e\a)],(34) 



where 



k' = Vl - P , 
e = arcsin^///i , 



(35) 
(36) 



where the index J physically means the phase-winding 
number. Henceforth, we call the condition p3l) the phase 



and Ao(e\Q;) is Heuman's lambda function whose def- 
inition is given in Eq. (|A20[) of Appendix A. Equa- 
tion ([M)) has a unique real solution k E [0,1] if and only 
if the soliton-train solution exists in the (7, f2)-parameter 
plane. Otherwise only the plane-wave solutions with an 
arbitrary phase-winding number exist. That is, when 
there is a real solution k of Eq. (|34p , the GP equation P7)) 
has the soliton solution in addition to the plane- wave so- 
lution, which is always a formal solution of the GP equa- 
tion. 

For repulsive interactions both the soliton and the 
plane-wave solutions are stable, as we demonstrate in 
Sec. IIVI However, for attractive interactions (7 < 0) the 
situation is different: above a certain magnitude of at- 
tractive interaction the plane- wave solutions become dy- 
namically unstable (25], |4^ , and hence the stability of the 
plane-wave solution does not hold for attractive interac- 
tions. The metastable bright soliton-train solutions can 
be obtained in a similar manner as shown in Appendix 
A. 

We solve the phase-boundary condition of Eq. for 
j = 1,2, and 3 solitons for repulsive interactions, and 
show the solutions in Figs. [^Ja) -[D^c), respectively. In 
the filled regions the soliton-train solution with j density 
notches coexists with plane wave solutions of arbitrary 
phase winding, and the unfilled region signifies that the 
soliton-train solution with j density notches does not ex- 
ist. The phase boundaries written in solid curves will be 
later identified as the parabolas given by Eq. (j56p from 
the Bogoliubov theory of Sec. IIV[ each boundary curve 
corresponds to different value of J in Eq. (|56p . The inte- 
ger value written next to each parabola in Fig. [Ha)-(c) 
represents different values of the phase-winding number 
J. 

We can analyze the phase boundaries in Fig. [2] as fol- 
lows. Two parabolic curves with two distinct phase wind- 
ings J and J' intersect at certain values of fi, which we 
denote by f2nodos(|>/ — J'\), satisfying \ J' — J\ = j. That 
is, a soliton state with j density notches and phase wind- 
ing J, and one with j density notches and phase winding 
J' become the same sn soliton with j nodes at f^nodcs, 
provided that the difference in J and J' is equal to j. 
Along lines ilnodos(|>/ — "^'D for arbitrary 7, the phase 
of the soliton slips by 2tt\J ~ J'\ and the wavefunction 
is given by the black soliton train with j density zeros, 
or nodes. For instance, for the triple-soliton-train case 
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j = 3 [Fig. Hfc)], in the filled region the soliton solu- 
tion ipj'=2 j=3' ^'^'^ j=3 exists in < 17 < 0.5, and 
0.5 < il < 1, respectively. There is thus a discontin- 
uous phase jump at r2nodes(| J — J'\) ~ 0.5, where the 
soliton solution is written as the black soliton train with 
three density zeros. The phase-jump lines correspond to 
17nodos(3) G {±0.5, ±1.5, . . .} for j = 3 for arbitrary 7. 
We note that the phase jump also occurs for plane waves 

at rinodcs- 

Typical behavior of the elliptic modulus fc, determined 
by Eq. ([M]) for j = 1, is shown in Fig.[5Jd). The solution 
k is zero on the phase boundary, rapidly grows once the 
parameters (7, f2) enter the soliton regime, and quickly 
approaches unity. This behavior is qualitatively the same 
for j > I soliton trains. 



this we study the energy and chemical potential of all 
stationary states. Derivatives of these quantities charac- 
terize the order of the metastable phase transition; our 
use of the term metastability refers to the fact that we 
consider excited as well as ground states. The energy per 
particle and chemical potential of the plane-wave state 
are given by 



,(pw) 



(n - Jf + 



2' 



(37) 
(38) 



C. Metastable Quantum Phase Transition 



We identify the metastable quantum phase transition 
in the mean field theory, as alluded to in Sec. H] To do 



For soliton solutions the energy per particle and chemical 
potential are 



,(st) 



(st) 



7 + 

37 
2 



[2,KE - (2 - fc2)A'2] + 
[3/^^; - (2 - k^)K'^] 



2K^ 



(39) 
(40) 



The elliptic modulus fc, which appears both explicitly in 
Eqs. ([39| -(j40 l) and in the complete elliptic integrals K = 
A'(fc) and E = E{k), contains f2, 7, J, and j implicitly, 
as described in Eq. p4p . 

The energy diagram shown in Fig. [3] summarizes the 
key result of mean-field theory. The figure plots all the 
stationary-state energies (choosing the zero of energy 
to be 7/2 for plotting convenience) for various phase- 
winding numbers J as a function of il, with 7 being 
fixed. Taking 7/2 as the zero of energy removes the triv- 
ial 7 dependence of E^f^\ and we can thus study an 
excess/shortage of solitonic energy from the plane- wave 
energy for a given strength of interaction. The plane- 
wave energies are trivial parabolas, and the different val- 
ues of J (written as integer values next to the parabolas 
in the figure) result in a discrete phase-winding number 
with respect to VI. Let us call the regime £ — 7/2 < 0.5^ 
the ground-state regime, 0.5^ < £ — 7/2 < 1^ the first- 
excited-state regime, 1^ < £ — 7/2 < 1.5^ the second- 
excited-state regime, and so on. 

We substitute the solution k of the phase boundary 
condition of Eq. dM]) for 7 = 1 into Eq. (UHl), and plot 
— 7/2 for various values of J and j in Fig. |31Ja). We 
find that soliton branches bifurcate from the plane- wave 
branches with the same winding number at the phase 



boundary, so that 

^(pw) _ g 



(st) 
J,3 ' 



5n£ 



(pw) 

J 



(41) 



This shows that a plane wave with phase winding J can 
be continuously deformed into a soliton with the same 
winding number without an energy jump. 

Furthermore, at the points rinodcs, where the gray soli- 
ton train becomes a black soliton train, the difference in 
adjacent phase-winding numbers of solitons equals the 
number of density notches j of the soliton train as shown 
before, and satisfies 



,(st) 



£ 



(st) 



(42) 



which, again, shows that there is no energy jump as- 
sociated with the self-induced phase slip caused by the 
soliton. Similar relations hold for the chemical potential: 



,,(pw) 



at the phase boundary. 



(st) 



(43) 



(44) 



at the phase-slip points ilnodcs- The first derivative of 
the chemical potential is a second-order cross derivative 
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of the energy, and is discontinuous: 



(45) 



at the phase boundary, where /i = d£/dN. Therefore 
the phase transitionis second order. All of this is true for 
both repulsive and attractive interactions. 

For repulsive interactions [Fig. [31 a)], the ground state 
is a plane wave located in the ground-state regime, with 
a ground-state phase- winding number of J = [f2 -I- 1/2J , 
where \x\ is the floor function which gives the integer 
closest to but not below x. We note that there is no 
soliton solution in the ground-state regime and the bifur- 
cations of soliton branches from the plane-wave branch 
exist only for excited metastable states. This is because 
the total energy increases upon the formation of soliton 
trains for repulsive interactions, as the density modu- 
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FIG. 3: (a) Energy diagram of stationary states for a fixed 
strength of repulsive (7 — 1) interaction. The zero of energy 
is taken as 7/2. The parabolas correspond to the energies of 
the plane wave, £!^P") _^/2, for various phase-winding num- 
bers. Other branches bifurcating from the parabolas denote 
the energies of soliton trains, -Ej***' — 7/2 for j — 1 (located 
in the first-excited-state regime), j = 2 (second-excited-state 
regime), and j = 3 (third-excited-state regime) dark solitons, 
respectively. Integers in the figures denote the phase-winding 
number J of soUtons, which is equivalent to that of the plane 
wave. The right panel enlarges the path C which connects 
the plane-wave and soliton branches in 0.3 < $1 < 0.7 (rough 
end points are indicated with the open circles), starting from 
the plane-wave branch, passing through the soliton branch 
between the phase boundaries (filled circles), and again goes 
back to the plane-wave branch, (b) Energy diagram of sta- 
tionary states for attractive (7 = —0.4) interaction. The j — I 
bright soliton is located in the ground-state regime, and j = 2, 
and i = 3 bright soliton train are located in the first, and the 
second excited regime, respectively. 



lation costs in both the kinetic and interaction energies. 
The lowest possible bifurcation thus starts from the first- 
excited-state regime, forming the upward swallowtail- 
shaped structure [4l[ where these two branches coexist, 
as can be seen in Fig.[3{a). The area of this swallowtail 
vanishes in the nonintcracting limit 7 = 0, and the area 
increases with increasing magnitude of interaction. 

Turning to attractive interactions [Fig. [3]Jb)], we ob- 
serve that the density modulation gains in interaction 
energy. The soliton branch thus has lower energy than 
the plane-wave branch, forming the downward swallow- 
tail structure, as seen in the figure. For small attrac- 
tive interaction of —0.5 < 7 < the ground-state 
branch is either the plane wave or a single bright soli- 
ton state with a ground-state phase- winding number of 
J = [fl + I/2J [42|- the two solution types are exclu- 
sive. For soliton solutions in 7 < —0.5, the gain in the 
interaction energy is always larger than the loss of the 
kinetic energy. The soliton branch thus separates away 
from the plane- wave branch, and the ground state always 
becomes a single soliton state {j = 1). In a similar man- 
ner, the j = 2 soliton branch, which is located in the first 
excited-state regime in Fig.[3Jb), separates away from the 
parabola for 7 < — 2, and the two lowest-energy states 
are given by j ~ 1 and j = 2 soliton branches. As shown 
in [42!| , the bright soliton solution approaches a Jacobi cn 
function at Onodos, and approaches a Jacobi dn function 
at the midpoint between two adjacent node lines finodos- 

The continuity of the first and second derivatives of 
energy with respect to a parameter are one way to iden- 
tify the order of a ground-state quantum phase transition 
driven by that parameter. We generalize this idea to the 
metastable two-parameter QPT by using the determi- 
nant of the Hessian matrix of a function f{Q,j), 



Det[i?(/)] = 



(46) 



From Eq. Bct[H [S'^^""'' )] = Det[i/(/i[f*^)] = for 

arbitrary (7, fl) for the plane-wave solutions. In order 
to calculate Eq. (|46)) for soliton solutions near the phase 
boundary, we numerically calculate the second deriva- 
tive of energy and chemical potential for each parameter 
(7, r2), and the first cross derivative with respect to both 
parameters. Figure |4] shows the structure of the deter- 
minant of the Hessian for the chemical potential /i, and 
energy £ in the (7, 51) plane. The result for the chem- 
ical potential diverges along the phase boundary, while 
the one for the energy is discontinuous along the phase 
boundary. The discontinuity for the latter increases as 
approaches flnodc, and it diverges at rinodo- 



D. Phase Winding and Unwinding 

Consider a condensate initially prepared in an excited 
metastable plane-wave state with a repulsive interac- 
tion. If it takes the continuous higher-energy path of 
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the swallowtail shown in Fig. ^a.) as is adiabatically 
changed, the metastable condensate undergoes an ener- 
getically smooth transition between distinct topological 
phases through a phase slip at rinodos- Figure [5] illus- 
trates the amplitude ^ p{9) and phase ip{9) along the 
path C (0.3 < O < 0.7) indicated in Fig. for a 

fixed repulsive interaction 7 = 0.6. At first the soliton 
branch does not exist, and the amplitude has a constant 
value of 1/V2tt. As the rotation fl increases, the ampli- 
tude starts to modulate at the phase boundary. The dips 
in the amplitude deepen as fl approaches r^nodcs = 0.5. 
At r^nodes the amplitude develops j = 3 nodes. At this 
point the wave function is described by the Jacobi sn 
function. Correspondingly, the gradient of the phase 
in the plane-wave regime is a constant value J, while 
the phase starts to wind in the soliton regime, keep- 
ing (p2.3id + Stt) — (p2.3{S) ~ 47r. For the sn soliton 
at D, = 0.5, the phase jumps by tt at the node posi- 
tions. For 57 > 0.5 the phase boundary condition be- 
comes ip-i^3{9 + 2n) — (^_i_3(6') = — 27r and the phase is 
gradually straightened, and gets unwound in the plane- 
wave regime. These changes occur in an energetically 
continuous manner. 

In a superfluid or metastable superflow, the angular 
momentum is quantized to be an integral multiple of N. 
In particular, the ground-state angular momentum L in 
a weakly repulsive ID superfluid ring system is quan- 
tized as an integral multiple of N at zero temperature, 
and there are discontinuous jumps between states having 
different values of the phase- winding number. In fact, 
this applies only to the ground state, and we show that 
the discontinuous jumps are replaced with the continu- 
ous crossover of angular momentum in metastable excited 
states. 

The average angular momentum per particle 



L 

N 







27r / a 



(47) 



of the plane- wave state is the integer J. In contrast, L/N 
of the soliton is nonintcger, as can be derived from our 
mean field formalism: 



L 



(st) 

J.J 



N 



S 



fgh 
2 



(48) 




FIG. 5: (Color online) Change in the amplitude |'i/'(^)| and 
phase ^p{9) in the third-excited-state regime for 7 = 0.6 for the 
higher-energy path C shown in Fig. [S] For ~ 0.3 the am- 
plitude has a constant value of \ip\ = l/\/27r and the gradient 
of the phase is equal to 2. As increases, a bifurcation in the 
energy occurs at f2 ~ 0.39 and both the amplitude and phase 
start to wind in the higher-energy soliton branch while keep- 
ing ip{27r) — f{0) = 47r. The phase-winding number changes 
from J = 2toJ = — latr2 = 0.5 and the amplitude has three 
nodes where the phase jumps by tt. For SI > 0.5 both the am- 
plitude and phase start to unwind with ip{2TT) — ip{0) = —2tt, 
and the winding disappears at SI ~ 0.61. The sequence of un- 
winding through a single soliton was described in our previous 
work [2411. 



As the second term becomes zero at the phase boundary, 
Eq. (|48p coincides with the average angular momentum 

of plane wave L^J^^ /N = J. Figure |6] plots the aver- 
age angular momentum L/N along the path C in Fig. [31 
for several strength of repulsive interactions. The phase 
winding and sign S are given by J = 2 and S = —1 in 
< f7 < 0.5, and J' = -1 and 5 = 1 in 0.5 < f7 < 1, 
respectively. We denote the critical angular frequency in 



each region by ri^rit' ^^"^ ^crif respectively. 

The critical frequencies are determined from linear per- 
turbation theory (Eq. ([5S[) in Sec.lIVIbelow) by using the 
values 7, 5, and the corresponding phase-winding num- 
ber in each region. As illustrated in Fig. [51 the angular 
momentum is smoothly connected aiVL = finodcs, and lin- 
early depends on VL in the soliton regime 
with a gradient of 



(2) 



crit < ^ < ^crit 



J' -J 



(2) 



(1) 



(49) 



The angular momentum is thus well fitted by a single line 



(a) Det [H(e)] 
6x10^ 




(b) Det [H(m.)] 



-1x10' 




7-(st) 

N 



= an + b 



in ni'X <n< nl^i. We find 



L 



(st) 



N 



J 



(50) 



jf7-(j + j')V'(i)' + i 



(51) 



FIG. 4: Evidence of a second-order metastable quantum 
phase transition: determinant of Hessian matrix for (a) en- 
ergy and (b) chemical potential. 



where J and J' denote the adjacent phase- winding num- 
bers that meet at 51,iodos, and b is determined by the 



condition afl 



(1) 



6 = J at the phase boundary. 
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Thus, we have three critical points in our primitive unit 
cell for fixed 7(> 0): there is a point at which the higher- 
energy soliton path begins, Part way through this 
path, a soliton train makes a transition from gray to 
black, at the critical point rinodcs, which is always either 
half-integer or integer and is independent of 7. Finally, 
the higher-energy soliton path ends at A simplified 
account of this sequence can be found in |2J] . 



IV. LINEAR STABILITY OF METASTABLE 
STATES 



Although all the stationary solutions of the GP equa- 
tion are listed in Sec. IIIIl they may or may not be sta- 
ble in response to perturbation. In this section we con- 
sider linear perturbation. We show linear stability for the 
two kinds of stationary solutions for repulsive interaction 
by studying fluctuations around the stationary states of 
the GP equation, and argue that the soliton branches in 
Fig. [3] can be realized in practice. 

A stationary solution ip{9) of the GP equation under 
a small perturbation S evolves in time as 



{m 



KWe''"*]}, (52) 



where (u„, w„), and A„ are given as eigenstates and eigen- 
values of the Bogoliubov-de Gennes equations (BdGE), 
respectively, and n S |Z| denotes the energy- index of 
excitations. Recalling the structure of solutions in the 
Bogoliubov formalism |43} , for each eigenvalue A„ with 



. Y =0.21 
. Y =0.42 

■ Y =0.63 

■ Y =0.84 



7=2 



J'=-l 



0.5 



FIG. 6: Change in the average angular momentum L/N along 
path C for a triple soliton train in the third excited-state 
regime, where the phase- winding number is J = 2 (—1) for 
< O < 0.5 (0.5 < f2 < 1). The average angular momentum 
takes the value equivalent to the phase-winding number J 
in the plane- wave regime, and linearly depends on Q, in the 
soliton regime. 



positive norm. 







de [|m„(0)|2 - |i;„(0)|2] = 1, 



(53) 



there is also an eigenvalue A„ = — A„ with negative norm. 
The BdGE predict an infinite set of such solutions. One 
exception to this rule can exist. This exception cor- 
responds to Namhu-Goldstone modes. For instance, a 
black soliton is at rest on the ring with respect to any 
background superflow in the rotating frame. The Gold- 
stone mode of the soliton corresponds to a center-of-mass 
translation of the soliton in this frame. Henceforth we 
consider only the Goldstone mode and eigenstates that 
satisfy (j53p , since the eigenvalues with negative norm do 
not have physical meaning. For the Goldstone mode the 
corresponding eigenvalue is zero, and the latter eigen- 
states that have positive norms can be real (i.e., positive 
or negative) or complex depending on the stability of the 
condensate mode. In general. An, £ C. 

Let us recall the excitations of a uniform superflow. 
When a plane-wave state with phase-winding number J 
is regarded as a condensate mode, fluctuations from that 
condensate mode are given by the eigensolutions of the 
BdGE with positive norm. 



A 



(J,pw) 



^/W(p + 2i)-2l{n- J), 



ui (X e 



i{j+i)e 



vi (X e 



-i(j-i)e 



(54) 
(55) 



where I e Z denotes the single-particle angular momen- 
tum of the excitation, which serves as a good quantum 
number since [/f (fi), L] = 0. 



For 7 > —0.5, all the eigenvalues A 



( J,pw) 



are real, as 

apparent from Eq. ([M]) . From Eq. (|54p we also find that 
several negative eigenvalues (associated with eigenstates 
of positive norm) appear when we take a metastable ex- 
cited state as a condensate mode. These negative eigen- 
values correspond to other plane-wave branches located 
in lower energy regimes than the input condensate mode 
itself. For the case of repulsive interactions, the number 
of negative eigenvalues thus coincides with the number 
of metastable states that are located in a lower energy 
regime than the metastable state under consideration. 
In Fig. [7] we show excitation energies with respect to the 
plane-wave state with J = 2 in the third excited-state 
regime for < r2 < 0.5. For 0.5 < < 1 the excita- 
tion energies from the plane- wave state with J = — 1 are 
symmetric with respect to O = 0.5. 

One of the negative eigenvalues changes its sign at a 
certain set of parameters (7ciit, f^crit)- This set of critical 
values is found by equating I to be Sj in Eq. (|54[) and 



imposing the conditionA 



(J,pw) 
l=Sj 



as. 




7crit 



(56) 



This equation has two real solutions 5 = ±1. These solu- 
tions determine the phase boundary at which the soliton 
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FIG. 7: Eigenvalues of the BdG equations for a fixed strength 
of repulsive interaction (7 = 1). The input condensate mode 
(zero excitation energy shown in the solid line) is taken as 
the higher-energy metastable state in the third excited-state 
regime of Fig. |3ja). Dashed and thick curves plot excitation 
energies from the plane wave and soliton train states, respec- 
tively. 



branch starts or ends its coexistence with the plane wave. 
The region 7 > 7crit is identical to requiring that the 
phase boundary condition Eq. (p4)) has a real solution. 

In Fig. [7] arc shown the Bogoliubov excitation energies 
from the gray soliton branch '!/'j=2 j=3 ^lo^ig half of the 
path C indicated in Fig. [3] The eigenvalues of the BdGE 
with the soliton condensate mode taken as the station- 
ary state in Eqs. (|52p are also real in the attractive 
case, so long as 7 > —0.5. The soliton-train solutions arc 
therefore linearly stable (they arc also nonlincarly sta- 
ble [4^ for 7 > 0, although we do not demonstrate that 
here, and, according to a quantum tunneling analysis, 
metastable for — 1 < 7 < with an exponentially long 
lifetime [i^, where the difference in the critical point 
is associated with the inclusion of the Fock term). The 
excitation energies from the soliton branch arc found to 
be close to those from the plane-wave branch. The no- 
table feature in the soliton regime is that there appears 
a Nambu-Goldstone mode, which is continuously con- 
nected with one of the negative eigenstates with I = Sj 
from the plane wave. This mode reflects the sponta- 
neous symmetry breaking of the soliton-train state. At 
the point rinodos, a degenerate pair of excitation branches 
emerges, where the phase jumps up or down by tt at each 
soliton in the soliton train. For i7 > i^nodos the excitation 
branches from the soliton train are also symmetric with 

respect to r^nodcs- 

We proceed to consider the attractive case in more de- 
tail. While all of the excitation energies from the plane- 
wave state A/ are real for any 7 > —0.5, some of the 
eigenvalues become complex for 7 < —0.5 as seen from 
Eq. ([M)) . Independent of f2, the excitation energies with 



I = ±1 (i.e., excitation modes of e''^''^'^-'^)^) become com- 
plex for 7 < —0.5, and those with I = ±2 (i.e., excitation 
modes of e**^-^^^^*) become complex for 7 < —2, indi- 
cating the modulational instability of these modes [46j . 
These complex modes indicate that in Fig.[3l^b) the j = 1 
soliton branch separates away down from the plane- wave 
branch in the ground-state regime for 7 < —0.5, and 
the soliton branch becomes the ground state while the 
plane-wave branch is no longer stable. Similarly, when 
the J = 2 soliton branch separates from the parabolic 
plane- wave branch, it becomes the first excited state; the 
plane-wave branch has complex A; and the plane wave 
and soliton cannot coexist, in contrast to the repulsive 
case. 

For the potential experimental verification of the con- 
tinuous crossover between the distinct topological states, 
one may use circular waveguides or toroidal traps to con- 
fine a weakly interacting atomic cloud. First one has 
hot atoms above the condensation critical temperature, 
subjecting to a rotating drive with a certain angular 
frequency. In order to set a system to a metastable 
uniform superflow, one quickly stops the rotation and 
then lowers the temperature to make the cloud condense. 
Then the angular frequency of the rotating drive should 
be changed adiabatically. In this process, microscopic 
roughness or small distortion of the trap is enough to 
create noise sufficient to break the translation symmetry 
of the condensate, making it take the higher-energy path 
of a metastable soliton state, without any artificial dis- 
torting the trap. Finally, one stops the adiabatic change 
in the frequency, results in a superfiow with a different 
winding number from the initial state. 

Supposing parameters approximately those realized in 
the experiment of Gupta et al. [Hi, using N = 10^ ^'^Rb 
atoms in a trap with a transverse frequency 27r x 50 
Hz of the circular waveguide with the radius i? = 1 
mm, the dimensionless coupling constant with the de- 
fault three-dimensional scattering length of rubidium is 
(7id = 27r X 5.5 x 10~^. This results in the mean-field 
interaction strength 7 ~ 3. Although this value is about 
five times larger than the effective strength of interaction 
where the Gross-Pitaevskii type of mean-field approxima- 
tion is quantitatively valid, the transitions between two 
topologically distinct states yet appear at ilcr whose or- 
der of magnitude is determined by Eq. (|56|) . The initial 
angular frequency of the rotating drive fi is arbitrary. 
Once the persistent current with arbitrary J, which is 
determined by the initial fi, is fixed by stopping the ro- 
tation, the number of solitons j is automatically deter- 
mined by Eq. ([55]) . 

Summarizing the results of our two theoretical meth- 
ods, mean field (GPE), and first-order quantum fluctua- 
tions (BdGE) for repulsive interactions, we have shown 
that the excitations from the plane wave in the j"^ 
excited-state regime have j thermodynamically unsta- 
ble modes. However, there is no modulational instability 
for repulsive interactions, and thus both the uniform su- 
perflow and gray-soliton states are stable in the (7, fi)- 
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plane. At the phase boundary, all energies of the sta- 
tionary solutions and the eigenvalues of the BdG equa- 
tions continuously connect to the excitations in the soli- 
ton regime without any energy discontinuity. In particu- 
lar, when one of the excitation energies from a plane- wave 
metastable state changes the sign from negative to posi- 
tive, a soliton branch with the same phase- winding num- 
ber appears, and the Nambu-Goldstone mode appears as 
a result. 



V. CONCLUSIONS 

We have studied metastable excited states of the one- 
dimensional Bose g function of interaction strength 
and rotation, showing the stability of superflow in a ro- 
tating ring trap. The study of such a system is part of 
the classic investigation of superfluidity [13] , a study to 
which we have added insight for metastable states. 

In the weakly interacting regime, all stationary states 
and their energy diagrams can be obtained in mean- 
field theory. Although it was previously known that the 
one-dimensional nonlinear Schrodinger equation has both 
plane- wave and soliton-train solutions on a ring (ssl . |40| , 
we have pointed out that the energy diagram is char- 
acterized by the smooth bifurcation of a soliton branch 
from a plane-wave branch in the rotating frame. This is 
the key to the continuous change in the topology of the 



condensate wavefunction 241, which can be character 



ized as a self-induced phase slip. It is possible to adjust 
the phase winding and unwinding through the phase slip 
via a soliton train with nodes, by rotating the ring trap 
adiabatically and/or changing the interactions. 

For repulsive/attractive interaction, the soliton branch 
has higher /lower energy than the plane- wave branch with 
the same phase-winding number. At the phase bound- 
ary where bifurcations occur, we showed that the excited 
stationary states undergo a second-order quantum phase 
transition. 

Going beyond mean-field theory, we used Bogoliubov 
theory to examine the linear stability of these stationary 
states, and found that both the plane wave and soliton 
train branches are linearly stable for repulsive interac- 
tions, and past a critical value of the interactions, linearly 
unstable for the attractive case. 

The continuous change between topologically distinct 
states we have presented in the work can be observed in 
experiments as follows. The crossover in the metastable 
state can be realized starting from hot atoms confined in 
a fast-rotating circular waveguide or toroidal trap. By 
stopping the rotation and lowering the temperature, one 
can obtain a metastable superflow of the condensate of 
uniform density. One should then change the angular fre- 
quency of the trap adiabatically. In principle, this trap 
must have a small deformation so that the higher-energy 
path of the swallowtail (see Fig. is selected, but in 
practice, no deformation of the trap need be forced on 
the system, because an infinitesimal perturbation, as is 



unavoidable in experiments, is sufficient. As the angu- 
lar frequency is increased further across the degeneracy 
point where the gray solitons develop nodes and become 
"black" , the phase- winding number changes over the self- 
induced phase slip, and eventually reaches again a super- 
fiow of uniform density with a different phase winding 
from the initial state. In this way one is able to observe 
the phase winding and unwinding without any abrupt 
energy discontinuity. 

We make the conjecture that the qualitative features 
of our study and the connection between quantum phase 
transitions and semiclassical bifurcations can be found 
in many quantum field theories. For instance, we expect 
that ID pscudo-spinor bosons on a ring display the same 
kinds of features, with the emergence of nonlinear ob- 
jects in the form of spin-textures signaling a metastable 
QPT in the same way that dark solitons did in our scalar 
theory. Such features may also appear in fermionic the- 
ories where the energy-gap function takes the place of 
our mean field. It is an open question as to the validity 
of our concept in higher dimensions, since vortices are 
fundamentally quantized, unlike solitons; however, the 
presence of boundaries may provide for the same kinds 
of features as we have described in ID, since a vortex nu- 
cleating on such a boundary can gradually approach the 
symmetry axis of a system and thereby increase the aver- 
age angular momentum continuously. However, the pos- 
sibly of a vortex lattice complicates the matter. Again, 
in spinor theories we expect vortex textures to take this 
role. Finally, we point out that in the present exam- 
ple both the semiclassical mean-field and the underlying 
many-body Hamiltonian were integrable. It would be 
intriguing to consider an example of metastable QPTs 
in a system for which one or both of these limits were 
non-integrable. 

This material is based upon work supported by the 
Sumitomo Foundation (RK), the National Science Foun- 
dation under Grant No. PHY-0547845 as part of the 
NSF-CAREER program (LDC), and a Grant-in- Aid for 
Scientific Research (Grant NO. 17071005) (MU). 



Appendix A: Stationary Solutions of the Nonlinear 
Schrodinger Equation 

In this appendix we provide a detailed derivation of 
soliton-train solutions of the GP equation for both repul- 
sive and attractive interactions. Substituting the general 
form of the solution ip{9) = y^e'-^C') into the GP 
equation, and equating the real and imaginary parts re- 
spectively, we get 



(Al) 



VP 



. (A2) 
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By integrating Eq. (|A1[) we have 

7r7p3 + Vp- -W^^ Mp'- 



(A3) 



The solution of this equation is given by 

= \ ; ^ ^ (A4) 

A/'(?7)Wdn2 fi^^^^^,fcj -7?/fc'2, 7<0 

where k"^ + k''^ = 1. The normahzation constant J\f is 
determined from d9p{9) = 1 as 



The depth 77 of the density notches is obtained from sub- 
stitution of Eq. (KM into Eq. (|XT|) : 



-20-if)Vg e [-1,0], 7>0, 
- 1 .9/[2(jfc'iC)2] e[0,l], 7<0. 



(A6) 



The integral constant W from Eq. (jA3p is given by 

where we defined functions /, g, h, and 5 for notational 
simplicity as 

/ = ±[Tr^j-2UKf + 2fKE], (A8) 
g = ^''^ + 2fKE, (A9) 
h = ±[TT^-f-2{jKy + 2fKE + 2{jkKy],{AlO) 

where the ± sign is for repulsive/attractive interactions, 
and 



S = sign(17 — J) 



+1, n> J, 



-1, n< J. 

From Eq. ()Aip we obtain the chemical potential 



(All) 



A* =^7+ - [3KE-{2-k')K']. (A12) 



I 

Note that there is no notational difference between the 
repulsive and attractive case in Eq. (jA12l) . By calculating 
the interaction energy per particle 

r2Ti 



Jo 



l-^fl] \3E^ -2(2- k^)KE + K\l-k^ 
2 87 \7r/ 



which is applicable to both repulsive and attractive cases. 

We next study Eq. (|A2p to obtain the phase pref- 
actor (p{0) and rewrite the phase boundary condition 
ip{9 + 27r) = ip{6) + 211 J. Equation (|X2|) can be read- 
ily integrated, giving 



W 

p 



(A15) 



By integrating this equation one more time, the phase 
part is obtained as 



(A16) 



where we used the definition of the elliptic integral of the 
third kind, 

n(e;u,fc) =y rfu[l -^sn^u]-^ (A17) 



Note that the parameter 



f 



(A18) 



is always positive (negative) for the repulsive (attractive) 
case. This difference in the sign is important to rewrite 
the phase boundary condition. 

Since the elliptic integral of the third kind becomes 
complete at 6* 27r as U{C,u,k) 2jU{^\a) [H, the 
phase boundary condition (p{0 + 2tt) — ip{d) + 2nJ is 
simplified for the repulsive and attractive cases as 

2 niVl - J)S 

_( 2{]k'Kf^2fl{gh) + v/27Vg + [l-Ao(e\a)], 



v/2^ + j^ [l-Ao(eV)], 
respectively, where 



arcsin ^/f/h, 
arcsin yjh/{k''^f), 7 < 



7 > 



(A19) 



one finds the expression for the energy per particle 



(A13) 



-J-J 



7 + 



[iKE - (2 - k^)K'^] 



3^ y^] [3£;^-2(2-fc2)ii:£;+A'2(l-fc2)] 



2K^ fj 



(A14) 



and a is given in terms of the elliptic modulus by 
a = arcsin(fc) [4^. The function Ao is called Heuman's 
lambda function and defined as 

Ao(e, k) = - [KEie, k'^) - [K - E)F{e, k'^)] , (A20) 

TT 

where F(u, k) and £'(u, k) are incomplete elliptic inte- 
grals of the first and second kinds, respectively. 
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